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EFFECT OF EDDY DIFFUSIVITY ON WIND-DRIVEN CURRENTS 
IN A TWO-LAYER STRATIFIED LAKE 
by Richard T. Gedney, W. Lick,* and Frank B. Molls 
Lewis Research Center 

SUMMARY 

The steady- state wind- driven circulation is numerically calculated in a rectangular 
stratified lake. The lake is composed of two layers (an upper, warmer layer termed 
the epilimnion and a lower, cooler layer termed the hypolimnion) having uniform but 
unequal densities and eddy diffusivities. The position of the interface between the two 
layers (the thermocline) and the three-dimensional velocities in both layers are calcu- 
lated using the shallow lake equations of P. We lander. Continuity of velocity and shear 
stress are maintained across the thermocline. 

Because the values of the eddy diffusivities in the layers are not well known, re- 
sults for hypolimnetic diffusivities from 1.05 to 16. 8 square centimeters per second and 
epilimnetic diffusivities from 16.8 to 67.2 square centimeters per second are given. 

In the examples calculated the length of the lake is of order 100 kilometers, the wind 
velocity is 5.2 meters per second, and the temperature difference between the epilimnion 
and hypolimnion is 18° C. The results show that, as the hypolimnetic diffusivity is 
increased, the thermocline tilt and hypolimnetic velocities increase. The effect of the 
other variables (wind stress, density, basin length, and mean thermocline depth) are 
easily discerned from the analysis . It is shown that the solution for a rectangular basin 
is not strongly dependent upon the length to width ratio of the basin. Therefore, most 
of the results are given for a square basin. 


INTRODUCTION 

As part of the present day analysis of the effects of pollution on large lakes, models 
are being developed which will attempt to predict the lake's chemical and biological 
♦Professor of Geophysics and Engineering, Case Western Reserve University, 
Cleveland, Ohio. Professor Lick's work was supported by the National Science Founda- 
tion. 




(a) Average temperature (b) Average temperature 
profile in Cayuga Lake profile in Cayuga Lake 
for period between Aug- for period between Sep- 
ust 14th and August tember 25th and October 
20th, 1968. (Data taken 1st, 1968. (Data taken 
from ref. 1. ) from ref. 1. ) 



Temperature, °F 

(c) Typical summer temperature profile in Lake Erie. (Taken 
from ref. 11. 1 

Figure 1. - Typical stratified lake temperature profiles. 

changes as a function of the pollution added to the system. Development of these models 
requires knowledge of the currents during the summer when the lake water has a very 
definite vertical temperature stratification. 

A typical thermally stratified lake, such as the one described in figure 1, can usu- 
ally be divided into three zones: the upper region called the epilimnion, the middle re- 
gion, where the temperature gradient is steepest, called the metalimnion, and the bot- 
tom region called the hypolimnion. The therm ocline is defined as the surface of 
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maximum rate of decrease in temperature. As is well known, the stable temperature 
gradient shown in figure 1 suppresses turbulence generated by the wind shear in the stir- 
face layers of the lake so that near the thermocline the turbulence intensity is greatly 
reduced from that encountered at the surface. This is well illustrated by the plot (fig. 2) 
of eddy diffusivity for heat as calculated by Sundaram, Esterbrook, Piech, and 
Rudinger (ref. 1) from the data which was used to construct the figure 1(a) temperature 
curve. In general, it has been found that the eddy diffusivities for heat and momentum 
under arbitrary thermal stratification conditions are a function of a stability parameter . 



0 100 200 300 

Eddy diffusivity of heat, v H , ftvday 


Figure 2. - Calculation of vertical eddy diffusivity in 
Cayuga Lake for period between August 14th and 
August 27th, 1968, (Data taken from ref. 1.) 


One of the more commonly used forms of the stability parameter is the Richardson num- 
ber Ri which is defined as 


Ri = a v g 


QT/9Z) 

(0U/0Z) 2 


where a y is the coefficient of volumetric expansion of water, g is the acceleration due 
to gravity, U is the horizontal velocity, Z is the vertical coordinate increasing upward, 
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and T is the temperature. Under conditions of neutral stability, Ri = 0, the eddy dif- 
fusivity of momentum and heat are generally assumed to be the same, that is, z^ M = 

The dependence of and v-^ on Ri has to date only been determined empirically. 
Because this empiricism has not been well established for v-^, the exact dependence of 
the turbulent Prandtl number v-^f on stability is a point of some controversy. It is 

generally believed for lake models (Hutchinson, ref. 2) that increases with in- 

creasing Ri. In any event, in the upper part of the epilimnion where the temperature 
gradient is small and as a result Ri is small, the ratio will be close to 1. 0. 

For the balance of the lake, v-^ should qualitatively have a similar profile as that 
shown in figure 2. 

The plots in figures 1 and 2 are based on average temperatures over a 1-week 
period during which time- dependent seiche currents existed. The internal waves asso- 
ciated with the seiche currents generate turbulence which diffuses the sharp temperature 
gradient in the neighborhood of the thermocline. During a period of steady winds, the 
temperature gradients in the hypolimnion and epilimnion and the thickness of the meta- 
limnion are generally much less than shown in the figures. 

The results just discussed for the temperature and eddy diffusivities suggest that 
the lake stratification may be modeled by considering the lake to be made up of two 
homogeneous layers each with different densities and eddy diffusivities with the inter- 
face between the two layers being located at the thermocline. At the interface, velocity 
and shear stress should be continuous. It is felt that this two- layer model which will 
be analyzed here will give a qualitative understanding of the dependence of stratified 
lake currents on the value of the governing parameters in the epilimnion and hypolimnion. 
A model in which the density and eddy diffusivities vary continuously with the vertical 
coordinate can be analyzed in the future when the dependence of v-^ on stability is 
better established. Since the two- layer model is most applicable under the condition of 
steady winds, it will only be solved assuming the wind stress is independent of time. 

Welander (ref. 3) and others have studied a two- layered model for the oceans which 
includes the Ekman dynamics assumption that the shear stress at the thermocline and 
the lake bottom are proportional to the geostrophic (inviscid) velocities in either the 
hypolimnion or epilimnion. In addition, at one point in their analyses the velocity in the 
hypo lim nion is assumed zero which allows the transport equation for the epilimnion to be 
uncoupled from the one for the hypolimnion. Others such as Hamblin (ref. 4) have used 
this approach in the Great Lakes. In fresh water lakes, the thickness of the friction 
layer generated by the wind is of the order of or greater than the average thickness of 
the epilimnion; and there is no geostrophic (inviscid) flow in the interior of the epilim- 
nion. Therefore Ekman dynamics cannot be used. In this analysis the shallow lake type 
equations originally derived by Welander (ref. 5) and shown by Gedney and Lick (ref. 6) 
to yield good quantitative results for Lake Erie during uniform temperature conditions 
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(no stratification) will be used. We also make no assumption pertaining to the hypolim- 
nion velocity magnitude and will determine the effect of the hypolimnion eddy diffusivity 
magnitude as well as the other governing parameters on the two- layer solution. 

Since the value of the eddy diffusivities are not well established, a wide range of 
values should be investigated. In order to accomplish this, a numerical solution of the 
complete two- layer equations was performed. This numerical solution is presented 
here. In a paper to be published an asymptotic expansion solution for the case when the 
hypolimnion eddy diffusivity is small will be given. 

Lee (ref.. 7) performed a numerical solution for a two- layer lake using the shallow 
lake type equations. Lee’s formulation and numerical procedure is different than the 
one used here and, because it allows the thermocline to intercept a variable bottom, is 
much more complex. His reported results, which are for one set of eddy diffusivities 
at wind shear stresses that may be too small, do not demonstrate the effect of eddy dif- 
fusivity. We consider here a parameter study over a wider range to demonstrate the 
effect of eddy diffusivity. In addition we are able to check our numerical results with 
analytical ones in the limit of small hypolimnion eddy diffusivities. 
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SYMBOLS 

coefficient defined by eq. (8) 
coefficient defined by eq. (9) 

wa^h + £) 

coefficient defined in appendix A 

coefficients defined in appendix A 

epilimnion friction depth 
hypolimnion friction depth 

coefficients defined in appendix A 

coefficients defined in appendix A 

Coriolis parameter 
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■Ill 


G ii 

coefficients defined in appendix B where i 

g 

acceleration of gravity 

H 

dimensional lake depth 

h 

nondimensional lake depth 

K, K h 

coefficients defined in appendix A 

L 

reference length of lake 

M 

number of x grid points 

M 1T 

Mfx + iM ly , epilimnion volume transport 

M 2T 

M2 x + iMgy, hypolimnion volume transport 

m 

unit normal vector to boundary 

m j, mg 

x and y components of m 

N 

number of y grid points 

_3_ 

JL + i_! 

3N 

3X 3Y 

_3_ 


3n 

3x 3y 

3 

J- - i-L 

3n* 

3x 3y 

P 

dimensional pressure 

Ri 

Richardson number 

T 

dimensional temperature 

U ref 

reference dimensional velocity 

u, U 

velocity in x direction 

v,V 

velocity in y direction 

w, W 

velocity in x direction 



y,Y \ 

Cartesian coordinates 

z,Z J 


a v 

coefficient of volumetric expansion 

“1 

2 ^ref /d l 


10 and j = 1, 2 
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a e 

t 

P 

T 


2 ^ref/ d 2 

length to width ratio in a rectangular basin 

Uj + ivj 

u 2 + iv 2 
mesh spacing 


? 

k t 

I 

P 

Ap 

^r 

T W 

U) 


dx* dy* 

nondimensional surface displacement 

Wr 

eddy diffusivity of heat 

eddy diffusivity of momentum 

nondimensional thermocline position 

density 

(P 2 " Pi)/Pi 

Pl/P 2 

wind shear stress t™ + ir^ 

x y 

(1 + 1)/2 


Subscripts: 

h indicated derivative with respect to h 

k indicates mesh point x location 

l indicates mesh point y location 

ref as in £ ref 

1 epilimnion 

2 hypolimnion 
Superscripts: 

(r) real part 

(i) imaginary part 

n iteration number 

( — ) indicates dimensional quantity 
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FORMULATION 


Basic Equations and Boundary Conditions 

In the present analysis, the lake is considered to be composed of two layers of dif- 
ferent density as shown in figure 3. In each layer, the basic approximations are that 
the water density is constant, the vertical eddy viscosity is independent of position, the 
pressure is hydrostatic (i.e. , it is determined by water column height), and the lateral 
friction and nonlinear acceleration terms can be neglected. The explanation for the den- 
sity and eddy diffusivity being considered constant in each layer has already been given 
in the INTRODUCTION. Gedney (refs. 8 and 9) has shown the other assumptions to be 
good approximations for the Great Lakes. 



Seiche currents which occur when the wind varies rapidly are very important in 
stratified lak es but there are periods when the winds are essentially steady. It is during 
the steady wind period that the temperature gradient at the thermocline is the sharpest 
and as a result the period when the lake is best approximated by two layers. The pro- 
perties of the eddy diffusivities are also known best for a steady wind. As a first step, 
the analysis performed here will assume steady winds. 

As is well known, the thermocline position in a stratified lake being acted upon by a 
steady wind will slowly sink. The rate of deepening has been measured in Cayuga Lake 
(see ref. 1) to be in the neighborhood of 15 to 30 centimeters per day. This rate is of 
such magnitude that the time derivative terms in the momentum equations can be ne- 
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glected and the problem is then "quasi- steady. " With the steady wind restriction plus 
the assumptions stated, the momentum equations and the continuity equation applicable 
to each layer shown in figure 3 are 


au. av. aw. 

— 2 +_i + 1 =o j = 1 , 2 

ax 3Y az 


- f c V ) 


+£ c u j = 


i sp j 

^ + v 

P] ax 


Mj 


a2 Uj 


az £ 


> j = 1, 2 


1 3P j 

+ v 

Pj av 


Mj 


a 2y j 

az 2 


a p . 

— 1 = -P^g j = 1, 2 

az •> 


(la) 


( 2 ) 


(3) 


As shown in figure 3 the Cartesian coordinate system used has X increasing eastward, 
Y northward, and Z vertically upward with the corresponding velocities being Uj, Vj, 
and Wj where j = 1 indicates the epilimnion and j = 2 indicates the hypolimnion. In 
these equations f is the Coriolis parameter, o- is the density, P. is the pressure, 
i^Mj is the vertical eddy diffusivity, and g is the acceleration due to gravity. Effects 
due to the Earth's curvature and to the variation in Coriolis force with position have 
been neglected since the scales of lakes are much less than the radius of the Earth. 

By integrating the hydrostatic equation (3) vertically from the surface of the lake, 

Z = £, to an arbitrary point in the epilimnion and from the thermocline, Z = £, to an 
arbitrary point in the hypolimnion, we obtain the equations for the pressure in the epi- 
limnion and hypolimnion as 

Pj(X, Y, Z) = pJ?(X, Y)\ - p lg Z + Pl gI(X, Y) 

P 2 (X, Y, Z) = pJjCX, Y)] - p 2 gZ + (p 2 - p 1 )gf(X, Y) + p lg ?(X, Y) 

In most cases the effects caused by the variation of P^ at the surface are small so that 
they can be neglected. Combining these pressure equations with the momentum equa- 
tions (2) results in 
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where 


‘Vmi 



X>l f c r l 

0N 


(4a) 


0 2 r 2 

P 2 U M2 T 
0Z^ 


- i P 2 f c r 2 = p lZ 


\ 0N 



(5a) 


Ap 


p 2~ p l 

P 1 


p r 


Pl 
p 9 


U 1 + 


iVi 


r 2 = U 2 + iV 2 
il = K + iK 

0N 0X 0Y 


1L = ^L + iiL 

0N 0X 0Y 


The aforementioned equations may be made nondimensional by introducing the fol- 
lowing variables: 
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4 = 


4 

^ref 


z 


Z 

^ref 


w =■ 


L W 
; ref ^ref 


u =- 


U 


U 

J 

ref 


v =- 


U 


V 

ref 


277Lt. 
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ref 


’ref 


Pi djg 


^ref 


£ref 
Ap ’ 
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nr. 


U 


_w 

ref 


ref 


p l d l f c 



Qf, = 


2 ^ref 


a r 


2 ^ref 



With these quantities the dimensional equations (la), (4a), and (5a) become 


3u t 0V-, 3w., 

— ± + — - + ± = 0 

0x 0y 0z 


(lb) 



0x 0y 


0w 2 
+ 

0Z 


= 0 


a i 



ir. 

2 ' 3n 


-i^-ir 2 =p r (ii + iS) 

a 2 dz 2 2 r U 3n/ 


(lc) 


(4b) 


(5b) 


where the continuity equation (la) has been written for both the hypolimnion and epilim- 
nion. 

This system of equations must now be solved subject to the following boundary con- 
ditions: 


<, 


w w , • w 1 

T = T x + lT y = — 


0T 




at z = 0 


a ^ 0z 


r 2 = w 2 = 0 at z - -h(x, y) 


r 1 - r 2 


,2 )■ at z=£(x,y) 

1 d 2 1 0r 2 


3z d 2 p r 0z 


( 6 ) 
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Wj = 0 at z = 0 


3 £ 0 £ ^ 

W 1 = U t — a + V, — 2- 

1 1 ax 1 ay 


> at z=|(x, y) 


'2 - 2 


9i ii 


3x 


9y_ 


( 7 ) 


where h is the depth of the lake. The and w 2 velocity boundary conditions at 
z = £(x, y) are such that there is no flow normal to the therm ocline. For strong enough 
winds the thermocline may intersect the surface of the water with an upwelling of cold 
bottom waters resulting. In this case, which will not be considered, there would be a 
velocity normal to the thermocline which would reflect the volume decrease in the hypo- 
limnion due to the upwelling. 

In this analysis all surface boundary conditions are applied at z = 0 (no wind posi- 
tion) instead of at z = £. For the cases to be calculated here it can readily be shown 
(see refs. 8 and 9 for error estimates for shallow Lake Erie) that evaluating the sur- 
face boundary conditions at z = 0 induces an error of a few percent which is only local 
in extent. 


Governing Equations for the Surface Displacement and the Thermocline 

Equations (4b) and (5b) subject to the boundary conditions (6) can be solved for the 
horizontal velocities in the hypolimnion and epilimnion as functions of the surface dis- 
placement £, the thermocline position i;, the lake depth h, and the wind stress t w . 

The resulting solutions are 

wai(^-z) cocMl+z) w ooa.z 

r . = A.,e 1 + A.e 1 +— e + 2i-2i £<z<0 (8) 

11 1 w 3n 

(9) 


<x>a 0 (ti+z) -va~(n.+z) 

: - A 2 e + 2ip r 


ii 

9n, 


-h < z < £ 


where 
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The Wj and Wg velocities can readily be determined, respectively, by combining 
equations (lb) and (8) and equations (lc) and (9). 

The velocity equations (8) and (9) and the continuity equations (lb) and (lc) can be 
integrated vertically to give 


M 1T - M lx + iM ly = 1° 01 l r 1 dz = D l^ + E 1 U + F 1 


3n ■ L 3n 


M 2T = M 2x + iM 2y 


^h 


2 r 2 d Z = D 2 T w + E 2 il + F 2 ii 

on on 


(10) 


( 11 ) 
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9M lx , 9M ly 


= 0 


3x 


3y 


3M 2x ^ dM 2y _ 


= 0 


dx 


3y 


(12a) 


(13a) 


where the coefficients Dj, Ep Fp Dg, Eg, and Fg are functions of the local depth 
h(x, y) and the thermocline position £(x, y), and are given in appendix A. The quantities 
Mjrp and Mgrp are respectively the volume transport in the epilimnion and hypolimnion. 
Equations (10) to (13a) can be combined to obtain two coupled equations for £ and 
Equations (12a) and (13a) are equivalent to 


Real 



= 0 


(12b) 


Real 



= 0 


(13b) 


where 


d _ _d_ _ 

3n* 3x 3y 


By substituting equations (10) and (11) into equations (12b) and (13b) we find the two equa- 
tions for £ and £ are 

F^V 2 ^ + E^V 2 £ + [e^— - + F^ — + F y — + F^ii + t w D^ - r w D^]ii 

1 J L 3x ^ 3y J h dx 3y ^ 3x x ^ 3x 


[4 


)!£ + E^)iS - F-P — + F^?P — + F^P-^- + r w D^ 


dx dy i h 


3x ^ dy ^ dy 


x 




M 

3y 


L(r) 3h E (i) 3 hi 3? _ f E (i) 3h _ E (r) 3 hi 3? _wf D (r) 3h D (i) 3 hi 

|>n p h 7; j'-'id'i; x L )h to i h ^J 

T y[ D jh , f- D ihf] = ° i - 1 - 1 
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where the superscripts on the coefficients indicate the real or imaginary parts and sub- 

j-t- 

scripts £ and h indicate the derivative of the j coefficient with respect to £(x, y) or 
h(x, y). The derivatives and so forth are given in appendix A. Equa- 

tions (14a) are solved subject to the boundary condition that the volume transport normal 
to the shore boundary is zero. Thus, if m^ and mg are the x and y components, 
respectively, of the outward unit normal m to the boundary, then for any closed body of 
water the boundary conditions are 

Mj x m 1 + M. y m 2 =0 j = 1, 2 (15a) 

Equations (14a) and (15a) are the boundary value problem for the solution to a com- 
pletely enclosed two- layer lake where the thermocline does not intersect the lake bottom 
or the surface of the lake. 


SOLUTION OF THE GOVERNING EQUATIONS 

In the previous section the mathematical model for the circulation in a two- layer 
lake was formulated. The model developed requires the solution of two coupled second- 
order nonlinear differential equations for the surface displacement £ and the thermo- 
cline position £. In a paper to be published, a partial analytical solution in the form of 
an asymptotic expansion will be given for the case when the vertical eddy diffusivity for 
momentum i^ M2 approaches zero. In this report we use finite differences to obtain a 
solution for equations (14a) subject to the boundary conditions (15a) for parameter values 
which do not permit use of the asymptotic solution. 






Figure 4. - Numerical grid for rectangular basin. 
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The point successive underre taxation iterative method was used to solve the govern- 
ing equations in a rectangular domain. The equations were linearized about the previous 
iteration by a Taylor series expansion. The linear finite difference equations at a grid 
point in the interior for the grid pattern shown in figure 4 are 



4? 


n+1 . ,.n+l 


k, l 


+ ^k, l - 1 + 




The G coefficients are listed in appendix B and are all of the form 


t n t n+l t n £ n+l -n *>n+l r n r n+i \ 

„ ^ ,^n «k+l, l " *k-l,Z *k, l+l “ *k, l-l ^k+l,l~ ^k-l,l *k, l+l ' ^k, l-l] 

G ij -°ij \^k,l’ — ” 


„n+l 


2A 


2A 


2A 


2A 


Here the previous and present iterations are indicated, respectively, by the superscripts 

2 

n and n+1. In these equations central differences with accuracy of order A were 
used for all derivatives except the first derivatives enclosed in the brackets [ For 

these derivatives either a forward or backward difference with accuracy of order A 
were used depending on which one increases the absolute value of the combined coeffi- 
cient for or Using forward or backward differences to increase the mag- 

nitude of the diagonal term in the coefficient matrix was necessary for the iteration 
scheme to be stable. 

For grid points along constant x and y boundaries, the finite difference equations 
for the equation (15a) boundary condition are 


E i 


(r) 


V 

+ f| p ) 
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Cl) 
IT'S , 

- 

[VI 

3x 
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L 9 yJ 

f/b J 
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8j 


0 = 1 , 2 
k - 1 
2 < l < 
x - 0 
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(15b) 
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L 3x Jb 


+ fM 


E i 


(i 


3? 


Lsyjf/b 


- F i 




J f/b 


+ G 7j*k|l - G 8j 


r \ = 1, 2 

k = M 

2<i<N-l 
x = 1 


(15c) 


r- 


E i 
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[3? 

+ F-{ r ) 

V 

+ E^ 

’s? 

+ F^ 

M 

L.3y. 

1 

f 

„3y. 

f 

_9x 

f/b 

_3xJ 


f/b 


- g 9^£;I = g ioj \ 


j = 1, 2 

1 = 1 

2 < k < M - 1 

ly = 0 


(15d) 
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(r) 


["3? 

+ F^ 


+ E^ 


+ f{^ 


i.3y_ 

] 

b 

.3y. 

1 

b 

_0X 

f/b 

.Sxj 


f/b 


+ G 9j^k|? " G 10j 


j = 1, 2 

1 = N 

2 < k < M - 1 
vj = 1 


(15e) 


where 


& 
.3x| f 


r n ~n+l 

_ ^k+1, l " *k, l 


i£ 

.3yJv 


^n+1 
_ ^k, l 


„n 


5 k- 1, l 


The Gg- to G 10j coefficients are given in appendix B and are functions of the n th 
iteration value of £ and £. For the derivatives tangential to the boundary in equations 
(15b) to (15e), either forward or backward difference was used depending on which one 
made the total coefficient of ^ and ^ maximum. Equations (15b), (15c), (15d), 
and (15e) were fised respectively for the corner points 1, 1; M, N; N, 1; and 1, M with 
the proper difference equation substituted for the first derivatives tangential to the bound- 
aries. 
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Equations (14b) and (15b) to (15e) constitute a set of equations for and at 

a mesh point. Every iteration, the equations were solved simultaneously for £ k + j and 
a t every grid point. The standard relaxation procedure was then applied with the 
relaxation coefficient usually taken to be 0. 5. In performing the calculations it was nec- 
essary to have at one point on the boundary a specific ? and £ value. In order to ac- 
complish this the difference between the new value at the boundary point and the desired 
value was either subtracted or added to the entire £ and £ field after each iteration. 

The iteration process was continued until the maximum ^ and £ relative error between 

_ ^ 

successive iterations at any point was less than 1 . 0x10 

It should be realized that the equations we are solving are highly nonlinear. The 
numerical procedure used to solve the equations (particularly the use of underrelaxation 
and the use of forward or backward differences to make the solution matrix diagonally 
dominate) was not known a priori but was developed empirically. Because the coeffi- 
cients in the equations that must be calculated every iteration are so complicated, each 
iteration takes a considerable length of time to perform. For a grid mesh of A = 0. 1 


(121 grid points over the lake region) a typical solution takes 20 to 40 minutes to obtain 
on an IBM 7090 system if the poor initial £ and C distributions of a constant are used 
to start the calculation. 


RESULTS 

The finite difference equations (14b) and (15b) to (15e) were solved for the position 
of the thermocline £ and the position of the lake surface t; . Using the solutions for £ 
and equations (8) and (9) determine the hypolimnion and epilimnion horizontal veloc- 
ities. 


Effect of Eddy Diffusivities 

In this section the effect of variation of the hypolimnion and epilimnion eddy diffu- 
sivities will be discussed. The results to be presented will all be for a constant depth 
square basin with the following parameter values being held constant: 
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Basin length, L, km 96. 6 

Basin depth, H, m 62. 6 

Epilimnion temperature, °C 22 

Hypolimnion temperature, °C 4 

A p 0. 00203 

p 0.9977969 

Wind velocity at 6 meters above water surface, m/sec 5.2 

t^, dynes/cm ^ -0.914 

Ty, dynes/cm^- 0.0 


A drag coefficient of 0. 00273 as proposed by Wilson (ref. 10) was used to determine the 
wind shear stress for the 5. 2-meter- per-second wind. The values listed previously are 
typical for the central basin of Lake Erie. 


Case 1 

In the first case to be discussed the eddy diffusivities used in the epilimnion and 
hypolimnion are respectively = 16. 8 square centimeters per second and = 

4. 2 square centimeters per second. These eddy diffusivities correspond to friction 
depths of dj = 18.2 meters and d 2 = 9. 1 meters. The wind drag coefficient and i/ M1 
value used here are the same as the ones used by Gedney (refs. 8 and 9) for current 
calculations in Lake Erie with uniform water temperature (no stratification). With 
these values for and the drag coefficient, the calculated velocities for Lake Erie 

agreed very well with measurements. The epilimnion in our model is considered to be 
of uniform temperature and for the cases calculated here will have an average thickness 
of approximately 20 meters. Since 20 meters is also near the average depth of Lake 
Erie, the value should be very similar to that used in the Lake Erie calculations. 

The value of v-y^ * n this ^ irs ^ cas e is somewhat arbitrary. Results for different values 
of v-yQ will b e reported in subsequent cases. 

Constant thermocline (£) contours when = 16.8 and Vy^ =4.2 are shown in 

figure 5(a). The thermocline assumes a shallow depth at the upwind end of the lake and 
a large depth at the downwind end. There is very little tilting of the thermocline in the 
cross wind direction. At a constant x station, the thermocline is deeper in the center 
of the lake than near the shores. The value of 9£/9n + 9£/9n throughout the basin is 
very close to zero. The thermocline therefore assumes a position so that the horizontal 
pressure gradi&nt 9 P 2 / 9 N in the hypolimnion is close to zero. With 9P 2 /9N being 
small, the velocities in the hypolimnion should also be small. 
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(c) Horizontal velocities at 3. 6-meter depth. (d) Horizontal velocities at 6. 7-meter depth. 

Figures. - Circulation in a two-layer basin of constant depth. Lake length, 96.6 kilometers; ■ 16.8 square centimeters per second; v M2 *4.2 
square centimeters per second. 
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(g) Horizontal velocities at 20-meter depth. 
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(f) Horizontal velocities at 15-meter depth. 
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(h! Horizontal velocities at 26-meter depth. 


Figure 5. - Continued. 




(i) Horizontal velocities at 32-meter depth. 
Figure 5. - Concluded. 


In figures 5(b) to (i), the horizontal velocities are given at the lake surface and at 
depths of 3. 6, 6.7, 10.0, 15.0, 20.0, 26.0, and 32. 0 meters. The dashed line included 
on some of the velocity plots is the intersection of the thermocline with the horizontal 
plane at that particular depth. The velocity patterns down to 10 meters are similar to 
those which occur in a homogeneous lake. The Coriolis force causes the deflection of 
the surface velocities to the right of the wind and a clockwise rotation of the current 
vector with depth. The effect of the epilimnion thickness being smaller at the upwind 
end of the lake causes the upwind velocities at the surface (fig. 5(b)) to be more in the 
direction of the wind. At depths closer to the thermocline, the epilimnion velocities 
rotate even more clockwise until they are in a southerly crosswind direction. As we go 
deeper to a region below the thermocline, the velocities become very small and aline 
themselves to the right of the wind. The velocity magnitudes at 32 meters are only 2 to 
3 percent of those at the surface. Note that the apparent increase in velocity between 
figures 5(f) and (g) only reflects a scale change. 
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(c) Horizontal velocities at 3. 6-meter depth. (d) Horizontal velocities at 10-meter depth. 

Figure 6. - Circulation in a two-layer basin of constant depth. Lake length, 96.6 kilometers; v^i - 67. 2 square centimeters per second; " 16.8 
square centimeters per second. 
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(e) Horizontal velocities at 15-meter depth. 
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(f) Horizontal velocities at 20-meter depth. 
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Figure 6. - Concluded. 







Case 2 


Results are given in figure 6 for the same parameters in ease 1 except and 

V M2 are Q. ua< ^ ru Pl e< ^ to be, respectively, 67.2 and 16.8 square centimeters per second. 
For comparison purposes the same thermocline depth at the upper righthand corner have 
been used in both cases 1 and 2. As shown in figure 6(a) the larger eddy viscosities pro- 
duce a greater crosswind thermocline tilt and maximum thermocline depth. As shown 
by figures 6(b) to (d) the velocities in the upper part of the epilimnion are approximately 
one- half that of the case 1 results. This is the same order of magnitude velocity change 
that occurs when the eddy diffusivity for a shallow (H/d^ < 1) homogeneous lake is quad- 
rupled. The significant result is that the velocities below the thermocline are larger 
than those in case 1 reflecting the increase in momentum transport across the thermo- 
cline because of the large eddy diffusivities. 


Case 3 

The effect of changing u M2 as v M1 is held constant is shown in figure 7 which is 
a plot of the thermocline depth along the y = 0 axis. With t^ M1 = 16. 8 square centi- 
meters per second, plots are shown for u M2 = 4.2 and 1. 05 square centimeters per 
second. We see, as u M2 — 0, the difference between the thermocline depth at x = 0 
and x = 1 becomes less. In addition, as v 2 - 0: (1) the. cross wind tilt of the thermo- 
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(a) Thermocline depth contours. 


<b) Horizontal velocities at 15-meter depth. 
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(c> Horizontal velocities at 20-meter depth. 


(d) Horizontal velocities at 26-meter depth. 


Figure 8. - Circulation in a two-layer basin of constant depth. Lake length, 96.6 kilometers; v M j - 14.8 square centimeters per second; “ 1-05 
square centimeters per second. 
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cline becomes smaller; (2) the thermocline depth distribution along x - 0 and x = 1 
becomes independent of y; and (3) the thermocline distribution as a function of x along 
y = 0 and y = 1 become the same. These characteristics are shown in the figure 8(a) 
thermocline contour plot for the case when Vy^ = 1*05 square centimeters per second. 

The velocities at 15, 20, and 26 meters are shown, respectively, on figures 8(b) to 
(d) when ia^ = 14.8 square centimeters per second and ia^ = 1-05 square centimeters 
per second. The velocities below the thermocline can be seen to be smaller than those 
in case 1 where = 14.8 and ia^ =4.2. 


Comparison With Asymptotic Limit 

As VyQ — 0 it is possible to construct an asymptotic expansion of the governing 
equations (10) and (11) in terms of the small parameter 


K 


T 


d 

d 


2 _ 1 _ 
l p r 


A zeroth- order solution for the thermocline was obtained and will be completely given 
in a future report. The zeroth- order thermocline position along the boundaries of a 
rectangular basin with t™ = 0 and constant < 0 is 


fr 


h, j = constant at X = L 




* H V (p2 ■ p i )g 


" 2L? X I- ,2 

x +Ut| at X = 0 


(16) 


2(X - L)t 
(P 2 - P x ) g 


+ U- 


at Y = 0, /3L 


where (3 is the width to length ratio of the basin. This solution predicts the zeroth- 
order thermocline position on the boundaries independent of /3, ^2’ 311(1 11 

-a! 2 /2(h+£) 

should be noted that the asymptotic solution is only valid when e « k . For 

V M2 — !- 05 square centimeters per second, this approximation is valid provided the 
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thermocline never approaches within a few meters of the lake bottom. The asymptotic 
solution interior to the boundaries does, however, depend on /3, Vy^, anc * the 

lake depth h. 

Included in figure 7 is a plot of equation (16). The numerical solutions approach 
the asymptotic results as v-y^ decreases. It should also be noted that the zeroth-order 
expansion gives the result that 


ii +i£ = o 

3n 3n 

- 2 

in the entire lake domain. The numerical solution gives 3£/3n + 3C/3n of order 10~ 
for ^2 = 1-05 square centimeters per second. This agreement between the zeroth- 
order asymptotic and numerical solution is good since the first-order term in the asymp- 
totic expansion would be of this magnitude. 

In figure 9 the thermocline positions for the zeroth-order asymptotic solution and 
the v 2 - 4.2 square centimeters per second numerical solution are given when the epi- 
limnion water volume are approximately equal in both cases. 



Figure 9. - Thermocline position at y = 0 for vjy^'O. 0 and 4. 2 square centimeters per second and = 16. 8 
square centimeters per second with equal epilimnion volume. 
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(c) Horizontal velocities at 10-meter depth. 
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(d) Horizontal velocities at 20-meter depth. 


Figure 10. - Circulation in a two-layer Dasin of constant depth with river flow. Lake length, 96.6 kilometers; u M1 - 14.8 square centimeters per second- . 4 2 
square centimeters per second. ' 
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(e) Horizontal velocities at 25-meter depth. 


Figure 10. - Concluded. 

Effect of Varying Other Parameters 

The effects of varying the hypolimnion and epilimnionjzddy diffusivities has been 
given in the previous sections. The effect of changes in r w , L, A p, and epilimnion water 
volume can produce large changes in therm oc line slope. The qualitative effects of these 
parameters are easily deduced from the asymptotic solution given in equation (16) and 
were confirmed by the numerical solution for the case when = 16. 8 square centi- 
meters per second and v ^ =4.2 square centimeters per second. 

Changing the lake depth produced no significant change in the thermocline position 
until the distance between the lake bottom and the thermocline is less than the hypo- 
limnion friction depth d 2 - When H + £ is made less than d 2 for a constant depth 
basin, the thermocline slope increases until the thermocline intersected the lake bottom. 
Including a bottom slope where H + £ was made small did not change this behavior. 

A plot of the thermocline position for the same parameters as case 1 shown in fig- 
ure 5 but including an east to west river flow of 5380 cubic meters per second is shown 
in figure 10. This river through- flow which is similar in magnitude to that in Lake Erie 
occurs in the epilimnion only. The main effect of the through- flow is to increase the 
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crosswind tilt of the thermocline. A through- flow of equal but opposite direction would 
produce an opposite crosswind tilt of the thermocline. Some velocities for the east to 
west through- flow are included in figure 10. The river flow effect on the velocities is 
quite small being most significant at the upwind end of the lake where the epilimnion 
thickness is the smallest. 


CONCLUSIONS 

The solutibn of the two- layer lake equations in a rectangular basin with a uniform 
5. 2-meter-per- second steady wind was solved for the case when the thermocline does 
not intercept the bottom or surface of the lake. The velocity and shear stress was made 
continuous across the thermocline. The density difference between the two layers was 
assumed to be 0.00203 gram per cubic centimeter. The following conclusions were 
determined from the results of the numerical solution: 

1. The major tilt of the thermocline occurs in the direction of the wind with the 
maximum thermocline depth being at the downwind end of the lake. 

2. In the Northern Hemisphere, a small crosswind tilt of the thermocline occurs 
with the depth of the thermocline to the right of the wind being less than that to the left. 

3. As the eddy diffusivity in the hypolimnion is decreased the thermocline tilt and 
hypolimnion velocities decrease. For the range of hypolimnion diffusivities investigated 

o 

(1. 05 to 16.8 cm /sec), the pressure gradient in the hypolimnion was always found to 
be small (of order 10"^). 

4. The shape of the thermocline was found to be weakly dependent on the length to 
width ratio and depth of the lake. As a result, most of the results are given for a 
square basin. 

5. As the hypolimnion eddy diffusivity is reduced, the numerical solution approaches 
an asymptotic solution which assumes the ratio of hypolimnion to epilimnion diffusivities 
to be much less than one. The zeroth- order asymptotic solution predicts the thermo- 
cline slope to be directly proportional to the wind shear stress and basin length and in- 
versely proportional to the density difference of the two layers and the local thermocline 
depth. 

6. Since the crosswind variation in the thermocline depth is small, the thermocline 
can be approximated by a cylindrical surface generated from a single curve. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, March 22, 1972, 

136-13. 
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APPENDIX A 


COEFFICIENTS FOR THE DIFFERENTIAL EQUATIONS 
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APPENDIX B 


COEFFICIENTS FOR THE DIFFERENCE EQUATIONS 


In this appendix the coefficients Dj, Ej, Fj, and their derivatives are all evaluated 
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APPENDIX B 


COEFFICIENTS FOR THE DIFFERENCE EQUATIONS 


In this appendix the coefficients Dj, Ej, Fj, and their derivatives are all evaluated 
7 and the local h value. 

9 L 


G. . = - tUdQ + fS? ^ + F,y 


J j4 


3x 




3y 


jl 


3x 


x~j£ y j£ jh 


9x 


jh 


ah 

3y 


G 0 - = E^-^ + E&^-^ + 2fJ^^ + T^D,y + + f£) — - FSi? — 


*2j - ^ 


3x 


It 


3y 


j£ 


3y 


x~j£ T J y"j| T "jh — - "jh 


3x 


G q • = E$) + e\ 1 ) + e[V — + Eiu — 

•3J J? 3x ^ 


3y 


J jh 


3x 


jh 


3y 


g„. = E^ii^ - eIP + e£) - e2? — 


4j 17 _ “K 


3x J h 3y 3x 


G,. = Fj?M| n + Ej?W 2 5 n + eK + f£) — + T x D iff ‘ t v D i% — 

^3 3£ j? 0 X 0y 3^5 g x x y 0 X 


fall *£? + E w °SH ♦ F g> ^ + T><‘) ♦ T?D<dl 
L ^ 3x 3y 3y x ^ 3y 

L(r) 3h F (r) 3 hi 3£ n L(r) 3h p (i) 3hl 3£ n 

P 1 **' MTy\-^ + [ F MTy ^^Jay 

i i + u 9h . E (i) ah i 

:Jl 7 n h| ^ i h «^J 8y 

wT D (r) 3h D (i) 3 hi w[" D (r) 3h _ D (i) 3hl 

x p h| * ih? y L )h5 * )h| 


E (i) 3h _ E (r ) 3h 

j»«^ » hi ^. 


+ T„ 
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G 6 j - F ]f 



+ U)a«! + E (i) 3il EH + r E Wii - E<*)Bli^ 

is 0 X 3s 0 y^ g x _ 3s 0 y 3s 3xJ 3y 


wDr) _3h D (i) 3hl T wr D (i) _3h _ D (r) dhl G t 

x p h 3x i h a y J y [> 3x i h aj 5j • 


where 


, 2 F (r) 0 2 tp( r ) 

E (r)_ 3E j E (r) _ 3 E j 

V~’ M 


ar 


ah a^ 


and 


„ t n £ n+l 

ai n %+i ,i ~ ^k- 1, i 


ax 


2A 


t n t n+l 

3| n l+l ~ % l-l 


3y 


2A 


..n t n+l fc n , £ n+l 4j; n 

y 2^n = %+M + %-M + Jk, l+l + _% l-l \ *?k, l 

A 2 


The ac n /3x, 3C n /9y, and V 2 ? n difference equations are the same as for 


G 7 - = r w D^ - AW + Eg)iC - + 

7j x y ]4 34 g x 34 3 y 34 0 X 34 3 y 




G 9j - M + t v D S«’ + >g )iC + B& )ie + F if — + 4 — 

9j X j? y ]? 34 0y 3s 0 X 3s 0 y 3s g x 


r _ AW - T w n( r ) + g £ n 
G 10j - T x D j T y D j +tj 9j 
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